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We present event-by-event viscous hydrodynamic calculations of the anisotropic flow coefficients V2 
to vs for heavy- ion collisions at the Relativistic Heavy-Ion Collider (RHIC). We study the dependence 
of different flow harmonics on shear viscosity and the morphology of the initial state. V3 and higher 
flow harmonics exhibit a particularly strong dependence on both the initial granularity and shear 
viscosity. We argue that a combined analysis of all available flow harmonics will allow to determine 
rj/s of the quark gluon plasma precisely. Presented results strongly hint at a value (t]/s)qgp < 2/4n 
at RHIC. Furthermore, we demonstrate the effect of shear viscosity on pseudo-rapidity spectra and 
the mean transverse momentum as a function of rapidity. 



I. INTRODUCTION 

Hydrodynamics is an indispensable and accurate tool 
for the description of the bulk behavior of a fluid. The 
equations of hydrodynamics are just the conservation 
laws, an additional equation of state and constitutive 
relationships for dissipative hydrodynamics. The idea 
that ideal hydrodynamics can describe the outcome of 
hadronic collisions has a long history. Applications to 
relativistic heavy-ion collisions have been carried out by 
many researchers (see [l|, for an extensive list of refer- 
ences). 

Fluctuating initial conditions for hydrodynamic sim- 
ulations of heavy-ion collisions have been argued to be 
very important for the exact determination of collective 
flow obscrvablcs and to describe features of multi-particle 
correlation measurements in heavy- ion collisions [3|-|22|. 
Real event-by-event hydrodynamic simulations have been 
performed and show modifications to spectra and flow 
from "single-shot" hydrodynamics with averaged initial 
conditions [Tt], [201422] . An important advantage of event- 
by-event hydrodynamic calculations is the possibility to 
consistently study all higher flow harmonics in the same 
simulation without the need for an artificial construc- 
tion of an initial eccentricity, triangularity, etc. This is 
particularly important for the computation of U4, which 
receives strong contributions from elliptical deformations 
of the initial state, and V5, which couples to triangular- 
ity from fluctuations and to the ellipticity of the collision 
geometry p2l |. Recent 3+1D viscous hydrodynamic sim- 
ulations have highlighted the role of fluctuating initial 
states also on electromagnetic observables f23j . 

Different v n depend differently on 77/s and the details 
of the initial condition, which is determined by the dy- 
namics and fluctuations of partons in the incoming nu- 
clear wave functions. In this work we present quantitative 
results on the dependence of V2 to v§ on both the shear 
viscosity to entropy density ratio 77/s and the granularity 
of the initial state, and compare to experimental data. 

This paper is organized as follows. In Section |n] we 
introduce the employed second order relativistic viscous 



hydrodynamic framework. The explicit form of the hy- 
perbolic equations in t-t\ s coordinates and the numerical 
implementation are presented in Section [Ml We discuss 
the initial condition for single events in Section IIVI and 
explain the freeze-out procedure in Section |Vj Finally, 
results are presented in Section IVI1 followed by conclu- 
sions and discussions in Section IVTT1 



II. VISCOUS HYDRODYNAMICS 

In [l[ we introduced the simulation MUSIC for ideal 
relativistic fluids and extended it in [20j | to include dissi- 
pative effects. 

In the ideal case, the evolution of the system, created 
in relativistic heavy-ion collisions, is described by the fol- 
lowing 5 conservation equations 



U t^- L id 



0. 
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(1) 
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where T-^ 



is the energy-momentum tensor and Jg is the 



net baryon current. These are usually rc-cxprcssed using 
the time-like flow 4-vector u M as 



T& v = (e + V)u»u v - VsT > 

J% = PBU» , 



(3) 
(4) 



where e is the energy density, V is the pressure, pb is 
the baryon density and g^ v = diag(l, —1,-1,-1) is the 
metric tensor. The equations are then closed by adding 
the equilibrium equation of state 



V = V{e,p B ) 



(5) 



as a local constraint on the variables. 

Historically, these equations have first been solved in 
a boost- invariant framework [24j , eliminating the longi- 
tudinal direction and assuming uniformity in the trans- 
verse direction. At RHIC the central plateau in rapid- 
ity extends over 4 units. Hence, as long as one is con- 
cerned only with the dynamics near the mid-rapidity re- 
gion, boost invariance should be a valid approximation 
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at RHIC, restricting the relevant spatial dimensions to 
the transverse plane. Much success has been achieved 
by these (2+l)D calculations (see references in [lj and 
[25|, [26| for thorough reviews). However, in order to an- 
alyze experimental data away from mid-rapidity, inclu- 
sion of the non-trivial longitudinal dynamics is essential 

The next step in improving relativistic hydrodynamic 
simulations of heavy- ion collisions is the inclusion of finite 
viscosities. In the first order, or Navier- Stokes formalism 
for viscous hydrodynamics, the stress-energy tensor is de- 
composed into 

TZ = T£ + S»» , (6) 

where T-^ is given by Eq. J3|) 

The viscous part of the stress energy tensor in the first- 
order approach is given by 

S" u = t) (v"u v + Vu" - ^A^W a u a ^j (7) 

where A^ v = g^ v — u^u v is the local 3-metric and 
V M = A^ LV d v is the local spatial derivative. Note that 
is transverse with respect to the flow velocity since 
A^ v u v = and u v u v = 1. Hence, w' 1 is also an eigen- 
vector of the whole stress-energy tensor with the same 
eigenvalue e. 77 is the shear viscosity of the medium, 
which we assume to be constant. 

This form of viscous hydrodynamics is conceptually 
simple. However, this Navier-Stokes form is known to 
introduce unphysical super- luminal signals [33l - l35| , lead- 
ing to numerical instabilities. The second-order Israel- 
Stewart formalism [36T - l38j avoids this super-luminal 
propagation, as does the more recent approach in [3|| . 

In this work, we use a variant of the Israel-Stewart 
formalism derived in [4(| , where the stress-energy tensor 
is decomposed as 

= T;7 + tt^ . (8) 
The evolution equations are 

d^V" = (9) 

and 

A^A^u'd^ = - — (tt^ - S^) - ^ v {d a u a ) . 

Ttt 3 

(10) 

When dealing with rapid longitudinal expansion, it 
is useful to transform these equations to the T-r) s - 
coordinate system, defined by 

t = t cosh r] s , 

z = Tsmhr) s . (11) 

We obtain the following hyperbolic equations with 
sources 

d a T$ = -d a ir ab + F b (12) 



and 

d a ( U a 7T Cd ) = -(l/T^)(TT Cd - S Cd ) + G Cd (13) 

where F b and G cd contain terms introduced by the coor- 
dinate change from t, z to r, g s as well as those introduced 
by the projections in Eq. (|10p . and tv is the relaxation 
time. 

Our approach to solve these hyperbolic equations relies 
on the Kurganov-Tadmor (KT) scheme [4l|, [42j , together 
with Heun's method to solve resulting ordinary differen- 
tial equations. 

III. IMPLEMENTATION 

As mentioned above, the most natural coordinate sys- 
tem for us is the r — rj s coordinate system defined by 
Eq. (fTTj) . In this coordinate system, the conservation 
equation d^J^ = becomes 

d T (T,J T ) + d v (rJ v ) + d rh J r <° = , (14) 

where 

J T = (cosh?7 s J° - sinh rj s J 3 ) , (15) 

J Va = (cosh rj s J 3 - sinh7? s J°) , (16) 

which is simply a Lorentz boost with the space-time ra- 
pidity rj s = tanh _1 (z/i). The index v and w in this sec- 
tion always refer to the transverse x, y coordinates which 
are not affected by the boost. Applying the same trans- 
formation to both indices in Eq. ([5]), one obtains 

<9 t (tT tt ) + d v (rT VT ) + d Va (T r,sT ) + T^" (17) 

+ 3 T (T7T TT ) + d v (T7T VT ) + d 7h {^ T ) + IT™' = , 

a T (rT r ^) +d v (TT vr >°) + d Va (T™>) +T T ^ (18) 
+ d T (Tir T ^ ) + d v (TTT v ^ ) + d Va (ir^ ) + ir Tr >° = , 

and 

d T {rT TV ) + d w (rT wv ) + d Vs (T r >" v ) (19) 
+ d T {TTT TV ) + d w {TTT wv ) + d Vs {n^ v ) = . 

These 5 equations, namely Eq. (|T4l for the net baryon 
current, and Eqs. (p~7l IT8l IT9|) for the energy and momen- 
tum, are solved along with Eqs. (flU)) for the viscous part 
of the stress-energy tensor, which in a more explicit way 
of writing read 

d c (u c n ab ) = - —u T ir ab + -A a, V'7r br - -A^wV" 
It t t 

- g cf n cb u a Du f - — - -ir ab d c ii c 
+ H ( ~-A a, >g b, >u T + -A a,1 g bT u T 

T n \ T T 

+g ac d c u b - u a Du b - ^A ab d c u c ^j 
+ (a b) , (20) 
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The relaxation time is set to 3r)/(e + V), in line with 
the approach in (43|. It was also shown in [5H that the 
dependence of obscrvables such as v-i on tv is negligible 
when including the term (4/3)7r^(<9 Q u Q ) in Eq. ([TBI) . 

To solve the equations we use the KT algorithm as ex- 
plained in [l| . In detail, we compute the first step within 
Heun's method for Eqs. {TU [HI [THJ HJl) , then the first step 
for Eqs. (j2T))) . proceed with the second step for Eqs. (|TH 
[T71 [T8l IT9|) using the evolved result for 7r a6 , and finally 
compute the second step for Eqs. (f20|) . This concludes 
the evolution of one time step. 

One major difference to the ideal hydrodynamic equa- 
tions solved in [l| is the appearance of time derivatives 
in the source terms of Eqs. (JTB CHI HH E0J>- These are 
handled with the first order approximation 



g( T n) = (g(Tn) - S(t„-i))/At , 



(21) 



in the first step of the Hcun method, and in the second 
step we use 



g( T n) = (ff*(T„+i) - 9(t„))/At, 



(22) 



where g*{r n +i) is the result from the first step. 

As in most Eulerian algorithms, ours also suffers from 
numerical instability when the density becomes small 
while the flow velocity becomes large. Fortunately this 
happens late in the evolution or at the very edge of the 
system. Regularizing such instability has no strong ef- 
fects on the observables we are interested in. Some ways 
of handling this are known (for instance see Ref . ji^ ) . 

In this study, when finite viscosity causes negative 
pressure in the cell, we revert to the previous value of 
ir^ v and reduce all components by 5%. This procedure 
stabilizes the calculations without introducing spurious 
effects. 



IV. INITIALIZATION AND EQUATION OF 
STATE 

To determine the initial energy density distribution 
for a single event, we employ the Monte- Carlo Glauber 
model using the method described in [46| to determine 
the initial distribution of wounded nucleons. Before the 
collision the density distribution of the two nuclei is de- 
scribed by a Woods-Saxon parametrization, which we 
sample to determine the positions of individual nucleons. 
The impact parameter is sampled from the distribution 



P{b)db = 2bdb/(bl 



(23) 



where 6 m i n and b max depend on the given centrality class. 
Given the sampled initial impact parameter the two nu- 
clei are superimposed. Two nucleons are assumed to col- 
lide if their relative transverse distance is less than 



D — y ONNpK , 



(24) 



where jv is the inelastic nucleon-nucleon cross-section, 
which at top RHIC energy of ^/s = 200AGcV is ctnn = 



42 mb. The energy density is taken to scale mostly with 
the wounded nucleon distribution and to 25% with the 
binary collision distribution. So, two distributions are 
generated, one where for every wounded nucleon a con- 
tribution to the energy density with Gaussian shape and 
width Co in both x and y is added, one where the same is 
done for every binary collision. These are then multiplied 
by 0.75 and 0.25, respectively, and added. 

In the rapidity direction, we assume the energy density 
to be constant on a central plateau and fall like half- 
Gaussians at large |?7 S | as described in [l[: 



e(rj s ) cx exp 



2a2 



fl(lls|-lat/2) 



(25) 



This procedure generates flux-tube like structures com- 
patible with measured long-range rapidity correlations 
|47h49| . The absolute normalization is determined by 
demanding that the obtained total multiplicity distribu- 
tion reproduces the experimental data. 

As equation of state we employ the parametrization 
"s95p-vl" from obtained from interpolating between 
lattice data and a hadron resonance gas. 



V. FREEZE-OUT 

We perform a Cooper-Frye freeze-out using 

E 1T = rl T ^ = * / /K*VK<2% , (26) 
d 4 p dyp T dp T d(p p 

where gi is the degeneracy of particle species i, and £ 
the freeze-out hyper-surface. In the ideal case the distri- 
bution function is given by 

(2iry exp((u' i p AI - Pi)/T FO ) ± 1 

(27) 

where [ii is the chemical potential for particle species 
i and Tfo is the freeze-out temperature. In the finite 
viscosity case we include viscous corrections to the dis- 
tribution function, f = fo + Sf, with 



5f = / (1 ± /o)pW«/3 



1 



2(e + -p)T 2 



(28) 



where it is the viscous correction introduced in Eq. ([8]). 
Note that the choice Sf ~ p 2 is not unique [5l| and a 
potential source of uncertainties particularly at the larger 
studied px- 

The algorithm used to determine the freeze-out surface 
E has been presented in It is very efficient in deter- 
mining the freeze-out surface of a system with fluctuating 
initial conditions. 



VI. ANALYSIS AND RESULTS 

While in standard hydrodynamic simulations with av- 
eraged initial conditions all odd flow coefficients vanish 
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by definition, fluctuations generate all flow harmonics as 
response to the initial geometry. We follow [l9|, where 
V3 is computed in a similar way to the standard event 
plane analysis for elliptic flow, and for each v n define an 
event plane through the angle 



1 (sm(n6)) 

ip n = - arctan — — 

n (cos(n</>)) 



(29) 



Note that here we do not weigh the average by pt as 
done in 0, [2(3] and [H]. This is because we focus on 
the low momentum part of v n in our analysis and hence 
want the event planes to be defined mainly using the 
low momentum part of the spectra. Differences between 
the different definitions are however small and lead to 
variations of v n on the order of one percent or less. 
The flow coefficients can be computed using 



v n = (cos(n(</> - ip n ))) . 



(30) 



When averaging over events we compute the root mean 
square y because we compare to data obtained with 
the event-plane method (see [53]). First, we present re- 
sults for particle spectra as functions of pt and ?y s . Pa- 
rameters were chosen in order to reproduce the experi- 
mental data for the spectra when including all resonances 
up to 2 GeV (and some higher lying resonances to be con- 
sistent with what is included in the employed equation 
of state). The used parameters can be found in Table 
|TJ Values for the maximal average energy density (in the 
center of the system) (£ ma x) are quoted for most cen- 
tral (0-5%) collisions. In addition, all parameter sets use 
?7flat = 4.8 and a v = 0.7. 





(To [fm] 


tq [fm/ c] 


<£max)[GeV/fm 3 ] 


Tfo [MeV] 





0.4 


0.4 


65.7 


150 


0.08 


0.4 


0.4 


57 


150 


0.16 


0.4 


0.4 


50 


150 


0.08 


0.2 


0.4 


57 


155 


0.08 


0.8 


0.4 


57 


145 



TABLE I. Parameter sets. 

Fig. Q] shows the transverse momentum spectra of pos- 
itive pions, kaons and protons compared to experimental 
data from PHENIX Q in 20-30% central events. In 
Fig. [5] we present a comparison of the computed charged 
particle spectrum for 77/s = 0.08 in 15-25% central colli- 
sions as a function of pseudo-rapidity rj p with experimen- 
tal data from PHOBOS 

With the employed parameters we achieve very good 
agreement when including all resonance decays. In gen- 
eral, it is computationally too expensive to include reso- 
nances up to 2 GeV for all calculations. Hence, for most 
presented results we restrict ourselves to including reso- 
nances up to the 0-meson only. This is a good approx- 
imation because pions dominate the flow of all charged 
hadrons and it is mainly the p- and u- mesons that mod- 
ify the pion distributions. Fig. [3] shows how the v n for 
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FIG. 1. (Color online) Positive pion transverse momentum 
spectrum for 20-30% central Au+Au collisions using rj/s — 
0.08 including resonances up to 2 GeV (solid) and up to the (j>- 
meson (dashed) compared to data from PHENIX [5J]. Results 
are averages over 10 single events. 
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FIG. 2. (Color online) Charged hadron spectrum for 15-25% 
central Au+Au collisions including resonances up to 2 GeV 
(solid, averaged over 10 events) and up to the (^-meson (dot- 
ted, averaged over 100 events) compared to data from PHO- 
BOS m. 



charged hadrons are affected by including different num- 
bers of resonances. Including more resonances reduces 
all v n , however, the quantitative effect is small. The re- 
duction is caused by the kinematics of resonance decays. 
When including more resonances, at midrapidity one in- 
creases contributions from resonances at greater rapidi- 
ties where the anisotropic flow is weaker. The influence 
of higher lying resonances on W3 appears to be larger than 
that on the other v n . 

Next, we verify that our results are not plagued by 
large discretization errors. Higher flow harmonics are 
sensitive to fine structures in the system and for the case 
of ideal hydrodynamics with smooth initial conditions it 
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FIG. 3. (Color online) Charged hadron vi to vs for rj/s = 0.08 
as a function of transverse momentum pr averaged over 10 
single events, including resonances up to the 0-meson (upper 
end of each band) and all resonances up to 2 GeV (lower end 
of each band). 



was shown in [l| that V4 is very sensitive to the lattice 
spacing if it is not chosen small enough. Fig. [4] shows 
v u (pt) for two different lattice spacings, our standard 
value of a = 0.115 fm and a larger a = 0.2 fm. Differences 
are within the statistical error bars from averaging over 
100 events each. 
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FIG. 4. (Color online) Charged hadron i>2 to v$ for ij/s = 
0.08 and 00 = 0.4 fm as a function of transverse momentum 
Pt averaged over 100 single events for lattice spacings a — 
0.115 fm (solid lines) and a — 0.2 fm (dashed lines). 



Because we are presenting the first (3+l)-dimensional 
rclativistic viscous hydrodynamic simulation, it is inter- 
esting to demonstrate the effect of shear viscosity on the 
longitudinal dynamics of the system, which in a (1+1)- 
dimensional simulation was studied in [EH . l57j . 

Fig-El shows the modification of charged hadron 
pseudo-rapidity spectra caused by the inclusion of shear 



viscosity. The shape of the initial energy density distri- 
bution in the longitudinal direction is the same for all 
curves, which were each averaged over 200 events. The 
normalization was adjusted to yield the same multiplicity 
at midrapidity in all cases. In the range 2 < 1 77^, | < 4 the 
pseudo-rapidity spectra are increased, for larger r\ v de- 
creased by the effect of shear viscosity. We checked that 
this effect is almost entirely due to the modified evolution 
when including shear viscosity. The viscous correction 
to the distribution functions Sf (|2"8)) only causes minor 
modifications. Additional information can be obtained 
by looking at the average transverse momentum (px) as 
a function of rapidity. We show in Fig. [6] that also (pt) 
increases at intermediate rapidities and decreases at the 
largest \y\. For this observable the effect of Sf is larger. 

The modification in the viscous case is caused by the 
following effect: Faster longitudinal fluid cells drag slower 
neighbors by the viscous shear coupling. Naturally, the 
inclusion of both transverse and longitudinal spatial di- 
mensions in the simulation is needed to allow for such 
coupling. Fast fluid elements are slowed down, slower 
ones sped up, decreasing the number of fluid cells with 
the largest rapidities, increasing the number at interme- 
diate rapidities. Further diffusion then distributes the 
momentum in all directions, explaining the increase of 
(pt) at intermediate rapidities. 
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FIG. 5. (Color online) Charged hadron spectrum for 20-30% 
central Au+Au collisions for different values of rj/s including 
resonances up to the 0-meson. 



In Fig. [7] we show the dependence of v„(pt) on the 
shear viscosity of the system. Results are averaged over 
200 single events each. For vi to V4 we compare to exper- 
imental data from the PHENIX collaboration obtained 
using the event plane method [5a |. One can clearly sec 
that the dependence of v n (pT) on rj/s increases with in- 
creasing n. To make this point more quantitative, we 
present the ratio of the p^-integrated v n from viscous 
calculations to v n from ideal calculations as a function of 
n in Fig. [51 While v 2 is suppressed by ~ 20% when using 
77/s = 0.16, W5 is suppressed by ~ 80%. Higher harmon- 
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FIG. 6. (Color online) Positive pion average pr as a function 
of rapidity y for 20-30% central Au+Au collisions from ideal 
and viscous (r]/s = 0.08) including resonances up to the (j>- 
meson. 



ics are substantially more affected by the system's shear 
viscosity than i>2 and hence are a much more sensitive 
probe of rj/s. This behavior is expected because diffu- 
sive processes smear out finer structures corresponding 
to higher n more efficiently than larger scale structures, 
and has been pointed out previously in [l8j . 

So far all results were obtained using initial conditions 
with a Gaussian width oo = 0.4 fm. We now study the 
effect of the initial state granularity on the flow harmon- 
ics by varying ao. Decreasing ao causes finer structures 
to appear and hence strengthens the effect of hot spots. 
This results in a hardening of the spectra as previously 
demonstrated in fl7| . Because we want to compare to ex- 
perimental data, we readjust the slopes to match the ex- 
perimental pT-spectra by modifying the freeze-out tem- 
perature (see Table HJ. 

Fig. [9] shows the dependence of w„(pt) on the value of 
do, which we vary from 0.2 fm to 0.8 fm. While is 
almost independent of ao , higher flow harmonics show a 
very strong dependence. In Fig.[TU]we present the depen- 
dence of the ^-integrated v n on the initial state granu- 
larity characterized by ao . 

Higher flow harmonics turn out to be a more sensi- 
tive probe of initial state granularity than i>2- While we 
are not yet attempting an exact extraction of rj/s using 
higher flow harmonics, our results give a first quantita- 
tive overview of the effects of both the initial state gran- 
ularity and rj/s on all higher flow harmonics up to v§. 
Comparing Figs. [7] and [9j we see that v^ipx) obtained 
from simulations using rj/s = 0.16 is about a factor of 2 
below the experimental result, and that decreasing ao by 
a factor of two does not increase it nearly as much. Note 
that Co = 0.2 fm is already a very small value given that 
we assign this width to a wounded nucleon. It is hence 
unlikely that a higher initial state granularity will be able 
to compensate for the large effect of the shear viscosity. 
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FIG. 7. (Color online) pr-differential i>2 to «5 from ideal hy- 
drodynamics (left), viscous hydrodynamics with rj/s = 0.08 
(middle), and rj/s — 0.16 (right). Results are averaged over 
200 events each. Experimental data from PHENIX [581 ]. 



Similar arguments hold for v^(pt). 

A detailed systematic analysis of different models for 
the initial state with a sophisticated description of fluc- 
tuations is needed to make more precise statements on 
the value of rj/s. It is however clear from the present 
analysis that the utilization of higher flow harmonics can 
constrain models for the initial state and values of trans- 
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FIG. 8. (Color online) Ratio of charged hadron flow harmo- 
nics in viscous simulations to the result from ideal hydrody- 
namics. Results are averages over 200 single events each. 



port coefficients of the quark-gluon plasma significantly. 
The analysis of only elliptic flow is not sufficient for this 
task, because it depends too weakly on both the initial 
state granularity and r\j s. 

We present V2 and v% as a function of pseudo-rapidity 
in Fig. [IT] The V2(rj p ) result from the simulation is flat- 
ter than the experimental data out to rj p ss 3 and then 
falls off more steeply. A modified shape of the initial 
energy density distribution in the ^-direction, the inclu- 
sion of finite baryon number, and inclusion of a rapidity 
dependence of the fluctuations will most likely improve 
the agreement. 

In Fig. [12] we show results of v u {pt) for different cen- 
tralities using rj/s = 0.08. Overall, all flow harmonics 
are reasonably well reproduced. Deviations from the ex- 
perimental data, especially of v^,(j>t) in the most central 
collisions indicate that our rather simplistic description 
of the initial state and its fluctuations is insufficient. Im- 
provements can be made by a systematic study with al- 
ternative models for the fluctuating initial state based 
on e.g. the color-glass-condensate effective theory (along 
the lines of 0). 

Finally, the higher flow harmonics integrated over a 
transverse momentum range 0.2 GeV < pr < 2 GeV 
are shown in Fig.[TJJ] as a function of ccntrality. V2 has 
the strongest dependence on the centrality because it is 
driven to a large part by the overall geometry. The odd 
harmonics are entirely due to fluctuations as we have 
discussed earlier, and hence do not show a strong depen- 
dence on the centrality of the collision. 



VII. SUMMARY AND CONCLUSIONS 

We have demonstrated that the analysis of higher flow 
harmonics within (3+l)-dimensional event-by-event vis- 
cous hydrodynamics has the potential to determine trans- 
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FIG. 9. (Color online) Differential V2 and «3 (upper panel) 
and «4 and i»5 (lower panel) in 20-30% central collisions using 
rj/s = 0.08 and varying ao- Results are averages over 100 
single events each (200 events for <tq = 0.4 fm). 



port coefficients of the QGP such as rj/s much more pre- 
cisely than the analysis of elliptic flow alone. We pre- 
sented in detail the framework of (3+l)-dimensional vis- 
cous relativistic hydrodynamics and introduced the con- 
cept of event-by-event simulations, which enable us to 
study quantities that are strongly influenced or even en- 
tirely due to fluctuations such as odd flow harmonics. 
Parameters of the hydrodynamic simulation were fixed 
to reproduce particle spectra both as a function of trans- 
verse momentum px and pseudo-rapidity rj p . The studied 
flow harmonics to v§ were found to depend increas- 
ingly strongly on the value of rj/s and also on the initial 
state granularity. This work does not attempt an exact 
extraction of rj/s of the QGP but our quantitative results 
hint at a value of rj/s not larger than 2/Att. The reason is 
the strong suppression of V3 to W5 by the shear viscosity. 
A higher granularity of the initial state counteracts this 
effect, but our results indicate that this increase is not 
large enough to account for rj/s > 2/4ir. We will report 
on a detailed analysis of higher flow harmonics at LHC 
energies and a comparison to the experimental data in a 
subsequent work. 
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FIG. 10. (Color online) Ratio of v„ with initial granular- 
ity characterized by the Gaussian width ao — 0.8 fm to 
the case with ao = 0.4 fm and ao = 0.8 fm, respectively. 
Results are for 20-30% central collisions using rj/s — 0.08. 
Averages are over 100 single events each. 



FIG. 11. (Color online) V2 and U3 as functions of pseudo- 
rapidity rj p compared to data from PHOBOS [SsJ. Aver- 
ages are over 100 single events each. 
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